% This function computes the conditional standard deviations of y using a local
% linearization at x_t
function out = simConStdY_linear(xSimFull,model,pruningOn)

numSim    = size(xSimFull,2);
selectxf  = ones(model.nx,1);
selectnx1 = zeros(model.nx,1);
selectnx1(1:model.nx1,1) = 1;
if model.orderApp == 3
    xSim      = xSimFull([selectxf;selectnx1;selectnx1]==1,:);
elseif model.orderApp == 2
    xSim      = xSimFull([selectxf;selectnx1]==1,:);
elseif model.orderApp == 1
    xSim      = xSimFull([selectxf]==1,:);
end
    
nx         = length(xSim(:,1));
eta2Full   = zeros(nx,nx);
eta2Full(1:size(model.eta,1),1:size(model.eta,1)) = model.eta*model.eta';
if pruningOn == 1
    [~,y_cu]= modelFit_pruned(xSim(:,1),model);
else
    [~,y_cu]= modelFit(xSim(:,1),model);
end
ny      = length(y_cu);
stdYsim = zeros(ny,numSim);
ySim    = zeros(ny,numSim);
epsValue = 1D-6;
for s=1:numSim
    x_cu = xSim(:,s);
    if pruningOn == 1
        [~,y_cu] = modelFit_pruned(x_cu,model);
    else
        [~,y_cu] = modelFit(x_cu,model);
    end
    ySim(:,s) = y_cu;
    gx = zeros(ny,nx);
    for i=1:nx
        x_p = x_cu;
        x_p(i,1) = x_p(i,1) + epsValue; 
        if pruningOn == 1
            [~,y_p] = modelFit_pruned(x_p,model);
        else
            [~,y_p] = modelFit(x_p,model);
        end

        x_m = x_cu;
        x_m(i,1) = x_m(i,1) - epsValue; 
        if pruningOn == 1
            [~,y_m] = modelFit_pruned(x_m,model);
        else
            [~,y_m] = modelFit(x_m,model);
        end
        gx(:,i) = (y_p-y_m)/(2*epsValue);
    end
    tmp = gx*eta2Full*gx'; %+ %model.Rw(1:ny,1:ny);
    stdYsim(:,s) = sqrt(diag(tmp));
end

% Expressing conditional stds in deviation from the steady state level
xfMean          = zeros(model.nx,1);
vec_xfxf        = (eye(model.nx^2)-kron(model.hx,model.hx))\reshape((model.eta*model.eta'),model.nx^2,1);
xsMean          = (eye(model.nx)-model.hx)\(model.HHxxtil*vec_xfxf+0.5*model.hss);
xrdMean         = zeros(model.nx,1);
x0_cu           = [xfMean;xsMean(selectnx1==1,:);xrdMean(selectnx1==1,:)];
gx              = zeros(ny,nx);
for i=1:nx
    x_p = x0_cu;
    x_p(i,1) = x_p(i,1) + epsValue;
    if pruningOn == 1
        [~,y_p] = modelFit_pruned(x_p,model);
    else
        [~,y_p] = modelFit(x_p,model);
    end
    
    x_m = x0_cu;
    x_m(i,1) = x_m(i,1) - epsValue;
    if pruningOn == 1
        [~,y_m] = modelFit_pruned(x_m,model);
    else
        [~,y_m] = modelFit(x_m,model);
    end
    gx(:,i) = (y_p-y_m)/(2*epsValue);
end
stdY0           = sqrt(diag(gx*eta2Full*gx'));
%stdYsimDevSS    = log(stdYsim./repmat(stdY0(:,1),1,numSim));

stdYsimDevSS    = (stdYsim-repmat(stdY0(:,1),1,numSim))./repmat(stdY0(:,1),1,numSim);


%% Output
out.ySim            = ySim;
out.conStdYDevSS    = stdYsimDevSS;
out.conStdY         = stdYsim;
out.std_conStdY     = std(stdYsim,0,2);
labelYields = {};
for i=1:length(model.maturities)
    labelYields{i} = ['$r^{(',num2str(model.maturities(1,i)),')}_t$'];
end
out.labely  = [model.labely,labelYields];

end